On the numerical integration of the Lorenz-96 model, with scalar additive noise, for benchmark twin experiments

May 25, 2021

Authors:

Colin Grudzien¹, Marc Bocquet² and Alberto Carrassi³

University of Nevada, Reno, Department of Mathematics and Statistics
CEREA, A joint laboratory École des Ponts Paris Tech and EDF R&D, Université Paris-Est, Champs-sur-Marne, France.
University of Reading, Department of Meteorology and NCEO

Outline

Twin experiments in toy models
Multiscale model reduction error and random dynamical systems
Stochastic integration schemes and bias in model statistics
Numerical benchmarks:

Strong and weak convergence
Ensemble forecast statistics
Filtering statistics

Conclusions and future work

Bayesian Data assimilation in the geosciences

Ensemble-based forecasting and Data assimilation (DA) are the prevailing modes of:

prediction; and
uncertainty quantification

in geophysical modeling.
DA combines simulations from a physics-based numerical model and real-world observations of a physical process.
Bayesian framework:

An ensemble-based forecast is a sampling procedure for the forecast-prior probability density.

The physics-based numerical model (and previous estimates of the state) encapsulate the Bayesian prior knowledge.

The output of DA is an estimate of a posterior probability density for the numerically simulated physical state, or some statistic.

Bias in the forecast-prior estimate will introduce bias into the posterior probability density.
Therefore, understanding and quantifying model bias has become a central discussion in ensemble-forecasting and DA.

Bias in twin experiments and numerical precision

Deterministic, biased-model:
- the numerical precision of the ensemble forecast can be substantially reduced without a major deterioration of the DA cycle's (relative) predictive performance¹.
- Model bias overwhelms errors introduced due to precision loss when the model-twin is resolved with low accuracy.
Random, unbiased-model:
- differences in statistics of model forecasts of stochastic dynamical systems are observed due to the discretization errors of low-order schemes.
- Frank & Gottwald develop an order 2.0 Taylor scheme to correct the bias in the drift in the Euler-Maruyama scheme in stochastically reduced model².
Our study:
- how and to what extent bias in numerical integration schemes for random dynamical systems also biases twin experiment and DA statistics?
- In what ways can numerical precision be targeted for unbiased prior estimates from ensemble forecasts?

1. Hatfield, S. et al. Choosing the optimal numerical precision for data assimilation in the presence of model error. Journal of Advances in Modeling Earth Systems, 10, 2177–2191, 2018.
2. Frank, J. and Gottwald, G. A. A Note on Statistical Consistency of Numerical Integrators for Multiscale Dynamics. Multiscale Modeling & Simulation, 16, 1017–1033, 2018.

Single layer Lorenz-96 model

Image of global atmospheric circulation.

Courtesy of: Kaidor via Wikimedia Commons (CC 3.0)

We consider a simplified 1-dimensional model of the atmosphere around a latitude circle.
We will suppose that this latitude circle can be discretized into \( n \) total longitude sectors;

each sector \( i \) is represented by a single state variable \( x_i \).

The classical Lorenz-96 model is given as, \( \frac{\mathrm{d}\mathbf{x}}{\mathrm{d} t} \triangleq \mathbf{f}(\mathbf{x}), \) where,

for each state component \( i\in\{1,\cdots,n\} \), \[ \begin{align}\large{f_i(\mathbf{x}) =-x_{i-2}x_{i-1} + x_{i-1}x_{i+1} - x_i + F}.\end{align} \]

The state variables \( x_i \) have periodic boundary conditions modulo \( n \), \( x_0=x_n \), \( x_{-1}=x_{n-1} \) and \( x_{n+1}=x_{1} \).
The term \( F \) in the Lorenz-96 system is the forcing parameter that injects energy to the model.
The Lorenz-96 model approximates geophysical behavior with:

external forcing and internal dissipation with the linear terms; and
advection and conservation of energy in the quadratic terms.

Two Layer Lorenz-96 Model

From: Wilks, D. Effects of stochastic parametrizations in the Lorenz'96 system. Quarterly Journal of the Royal Meteorological Society 131.606 (2005): 389-407.

The two layer Lorenz-96 model simulates coupled, ocean-atmosphere dynamics.
This is a common model to study stochastic model reduction³;

the effects of fast-variable dynamics on the slow variables are parameterized with a stochastic process.

This is justified mathematically due to the Central Limit Theorem:

in the asymptotic separation of the time scales for the fast and slow variables,
the effect of the fast variables will reduce to additive, Gaussian noise⁴.

For finite separation of scales, the white-in-time model error assumption is no longer valid;
- non-Markovian memory terms should be included⁵.

3. Vissio, G. and Lucarini, V.: A proof of concept for scale-adaptive parametrizations: the case of the Lorenz’96 model, Quarterly Journal of the Royal Meteorological Society, 144, 63–75, 2018.
4. Gottwald, G. et al. Stochastic climate theory. Nonlinear and Stochastic Climate Dynamics. 209–240. 2015. Cambridge University Press.
5. Demaeyer, J. and Vannitsem, S.: Stochastic Parameterization of Subgrid-Scale Processes: A Review of Recent Physically Based Approaches, in: Advances in Nonlinear Geosciences, pp. 55–85, Springer, 2018.

L96-s Model

Define the L96-s model as follows, \[ \begin{align} \frac{\mathrm{d} \mathbf{x}}{\mathrm{d} t} \triangleq \mathbf{f}(\mathbf{x}) + s(t)\mathbf{I}_{n}\mathbf{W}(t), \end{align} \] where

\( \mathbf{f} \) is defined as in the single layer Lorenz-96 model
\( \mathbf{I}_n \) is the \( n\times n \) identity matrix,
\( \mathbf{W}(t) \) is an \( n \)-dimensional Wiener process; and
\( s(t):\mathbb{R}\rightarrow \mathbb{R} \) is a measurable function of (possibly) time-varying diffusion coefficients.

Both the truth-twin and model-twin are simulated with the L96-s model.

Numerical model is unbiased in representing the true physics, but where the physics are intrinsically random.

Represents an idealized two layer Lorenz-96 model in which the separation of the time scales of the atmosphere and ocean is taken to infinity.

This is a perfect-random model model assumption.

We distinguish effects of numerical bias from bias in the random physical process model.

Strong versus weak convergence

We study numerical integration schemes which converge in the strong sense;

strong convergence measures the expected path-discretization errors over all possible Wiener processes starting at an initial condition.
This is the analogue of deterministic path-discretization errors.

If \( \mathbf{x}_\mathrm{SP} \) is an exact sample path, evolving with respect to a particular Wiener process;
and \( \mathbf{x} \) is an approximation of this path, discretized at a maximum step size of \( \Delta \);
loosely, we say the approximate \( \mathbf{x} \) converges strongly to \( \mathbf{x}_\mathrm{SP} \) at order \( \gamma \) if:

there exists a \( \Delta_0 \) and a constant \( C>0 \) such that for every \( \Delta < \Delta_0 \)

then the expected path-discretization error is bounded as \[ \begin{align} \mathbb{E}\left[\left\Vert \mathbf{x}(T) - \mathbf{x}_\mathrm{SP}(T)\right\Vert\right] \leq C \Delta^\gamma \end{align} \] where these are the states, evolved from time \( 0 \) to time \( T \).

Strong versus weak convergence – continued

Weak convergence measures the error in representing some statistic of the forward distribution,
- given the evolution an initial point Dirac-delta distribution over all possible realizations of the Wiener process.
If \( g \) is a \( 2(\gamma +1) \) continuously differentiable function of at most polynomial growth;

we say that \( \mathbf{x} \) converges weakly to \( \mathbf{x}_\mathrm{SP} \) at order \( \gamma \) if for all \( \Delta< \Delta_0 \), \[ \begin{align} \left\Vert \mathbb{E}\left[ g(\mathbf{x}(T)) - g(\mathbf{x}_\mathrm{SP}(T))\right] \right\Vert \leq C \Delta^\gamma, \end{align} \] for any such statistic \( g \).

Loosely, strong convergence measures the mean of the path-discretization errors, while weak convergence can measure the error in representing the mean over all paths.

Strong versus weak convergence in twin experiments

Note: integration schemes that converge in the strong sense also converge in the weak sense,

however, weak schemes aren’t guaranteed to converge in the strong sense.

For twin experiments, there is an important distinction between these modes of convergence:

The truth-twin should be generated by a simulation which is precise in the strong sense,

here we assume we have observations of a path that is consistent with the governing dynamics ⁶.

An ensemble forecast representation of the prior needs to converge in the weak sense alone;
indeed, the forecast represents the sampling of the prior density, and we do not need to ensure the precision of each path solution if the density is accurate.

Therefore, the model-twin should be evaluated in terms of the precision in the weak sense.

6.Hansen, J. A. and Penland, C. Efficient approximate techniques for integrating stochastic differential equations. Monthly weather review. 134, 3006–3014, 2006.

Integration schemes

We study three commonly used numerical integration schemes for stochastic differential equations (SDEs):

Euler-Maruyama – a simple extension of the deterministic Euler scheme, order 1.0 strong convergence in L96-s.
4-stage Runge-Kutta – a simple extension of the determinstic 4-stage Runge-Kutta scheme⁷, order 1.0 strong convergence in L96-s
Second order Taylor – an integration scheme based on the Taylor-Stratonovich expansion⁸, order 2.0 strong convergence in L96-s.

The derivation of the order 2.0 strong Taylor scheme for the Lorenz-96 model is nontrivial and has not appeared earlier in the literature to the authors' knowledge.

Because L96-s model has:

constant or vanishing second derivatives in the deterministic component;
periodic boundary conditions; and
additive scalar noise;

we can compute this scheme efficiently, with complexity growing linearly in the system dimension \( n \).

7. Rüemelin, W. Numerical treatment of stochastic differential equations, SIAM Journal on Numerical Analysis, 19, 604–613, 1982.
8. Kloeden, P. and Platen, E. Numerical Solution of Stochastic Differential Equations, Stochastic Modelling and Applied Probability, Springer Berlin Heidelberg, page 359. 2013.

Strong convergence benchmarks

Plot of strong convergence discretization error versus step size.

Grudzien, C., Bocquet, M., & Carrassi, A. (2020). On the numerical integration of the Lorenz-96 model, with scalar additive noise, for benchmark twin experiments. Geoscientific Model Development, 13(4), 1903-1924.

Estimated strong convergence discretization error in the vertical axis versus the time-discretization step size;

log-log base 10 scale.

Point estimates — average estimated discretization error over 500 independent ensembles of simulations.

Each ensemble — evolves initial condition with respect to 100 independent Wiener processes.
We simulate:
1. finely discretized path solution \( \mathbf{x}_\mathrm{SP} \) with error on \( \mathcal{O}\left(10^{-7}\right) \), and
2. a coarser approximation by one of the tested methods.

Strong convergence benchmarks – continued

Strong convergence \[ \begin{align} \mathbb{E}\left[\left\Vert \mathbf{x}(T) - \mathbf{x}_\mathrm{SP}(T)\right\Vert\right] \end{align} \] estimated as:
1. the mean difference of the finely discretized solution, versus
2. the coarse discretization over all Wiener processes in the ensemble.
Ensemble-mean estimates are Gaussian random variables distributed around the true expectation.
Slope (order of convergence) estimated with weighted least squares;

weights proportional to the inverse standard deviation of the ensemble realizations.

Strong convergence benchmarks – continued

Diffusion level \( s \) fixed for each panel.
All the orders of convergence are empirically verified.

However, the constant \( C \) in the bound \( C\Delta^\gamma \) has a large impact on strong discretization error.

Low diffusion regime:

order 1.0 Runge-Kutta scheme has discretization error comparable to order 2.0 Taylor scheme.

Performance of Runge-Kutta scheme varies between low and high diffusion.

Order of convergence is same in all diffusion regimes, the difference in performance is reflected in the constant \( C \).

Weak convergence

Weak convergence \[ \begin{align} \left\Vert \mathbb{E}\left[ \mathbf{x}(T) - \mathbf{x}_\mathrm{SP}(T) \right]\right\Vert \end{align} \] estimated as:
1. difference of means of the finely discretized solution, versus
2. the coarse discretization over all Wiener processes in the ensemble.
Effect of constant is more pronounced in weak convergence.
Low diffusion:

1.0 Runge-Kutta outperforms 2.0 Taylor.

Weak discretization error varies on order of magnitude between low and high diffusion.
High strong precision across all diffusion makes Taylor a choice for generating the truth-twin.
Runge-Kutta can generate ensemble forecast in the model-twin for weak accuracy.

Ensemble forecast statistics

We use the Taylor scheme as a benchmark in the following experiments.

Taylor is consistent across diffusion levels.

We study:

how ensemble forecast statistics of,
1. Euler-Maruyama, and
2. Runge-Kutta
differ from Taylor.

Ensemble forecast statistics – fine discretization

Plot ensemble forecast bias versus time.

500 initial conditions of L96-s model:

each initial condition is forecasted with an ensemble of 100 independent Wiener processes.

We compute the empirical, ensemble-estimated mean and ensemble-estimated spread of the forward distribution.
This is performed over each of the integration schemes:

Taylor — step size \( \Delta=10^{-3} \)
Euler-Maruyama — step size \( \Delta=10^{-3} \)
Runge-Kutta — step size \( \Delta=10^{-3} \)

Top panels — RMSD of ensemble means for:
1. Euler-Maruyama and
2. Runge-Kutta
versus Taylor.
Bottom panels — ratio of the ensemble spreads for:
1. Euler-Maruyama and
2. Runge-Kutta
versus Taylor.

Ensemble forecast statistics – fine discretization continued

Plot of ensemble forecast bias versus time.

Summary statistics versus forecast time — over 500 initial conditions:

median is plotted as a solid line;
inner 80 percentile is plotted shaded;
min/ max values are plotted in dashed.

Runge-Kutta

ensemble statistics have almost no difference from Taylor up to \( T\approx 3 \).

Euler-Maruyama

ensemble mean — rapid divergence
ensemble spread — divergence and presence of bias :

median spread ratio above 1.0 asymptotically.

Ensemble forecast statistics – coarse discretization

Coarse ensemble versus benchmark:

Taylor — step size \( \Delta=10^{-3} \)
Euler-Maruyama — step size \( \Delta=10^{-2} \)
Runge-Kutta — step size \( \Delta=10^{-2} \)

Runge-Kutta

ensemble mean — faster onset of divergence.
ensemble spread — increased variance of ratio.

Euler-Maruyama

ensemble mean — extreme divergence, especially in low diffusion
ensemble spread — divergence and presence of bias :

MINIMUM spread ratio above 1.0 asymptotically.

Runge-Kutta does not have extreme divergence of trajectories, and the spread remains unbiased asymptotically.

Summary of bias in ensemble forecast

Runge-Kutta is robust in generating ensemble forecast statistics.

Runge-Kutta is very precise in the weak sense.

Coarse time step \( \Delta=10^{-2} \)

unbiased spread compared with Taylor using a step size of \( \Delta=10^{-3} \).
Divergence of ensemble means matches the asymptotic divergence with time-step of \( \Delta=10^{-3} \).
Low precision numerics inreases variance of sample forecast statistics, but remains unbiased.

Euler-Maruyama

systematic biases in ensemble forecast statistics.
Asymptotic ensemble spread is systematically, artifically inflated.
Extreme divergence of trajectories;
- divergence over an order of magnitude difference compared with the Runge-Kutta scheme.
Indicates Euler-Maruyama biasing ensemble forecast statistics,

doesn’t yet demonstrate the effect of bias on DA twin experiment.
This is demonstrated in the following…

Taylor benchmark configuration – filter statistics

Plot of benchmark Taylor configuration filtering statistics.

Benchmark configuration:

truth twin — Taylor step \( \Delta=10^{-3} \)
model twin — Taylor step \( \Delta=10^{-3} \).

Ensemble Kalman filter (EnKF) — RMSE and spread.
Vertical axis: diffusion level \( s \).
Horizontal axis: variance of the error in the observations given to the ensemble Kalman filter.
\( N=10^2 \) samples
State dimension \( n=10 \);
Fully observed

Convergence of filtering statistics without inflation / localization.
RMSE and spread are computed as the asymptotic average over \( 2.5\times 10^{4} \) analysis cycles.
Analysis RMSE is comparable to the analysis spread;

lower than the standard deviation of the error in the observations.

Comparision:

Vary the ensemble integration between Runge-Kutta and Euler-Maruyama, with a step size between \( \Delta\in\left\{10^{-3},10^{-2}\right\} \).

Filter benchmark – Runge-Kutta model-twin, coarse discretization, Taylor truth-twin, fine discretization

Plot of coarse Runge-Kutta filter simulation versus Taylor benchmark.

Generate ensemble with Runge-Kutta
Compare:

difference in RMSE with benchmark;
ratio in spread with benchmark

Runge-Kutta step \( 10^{-3} \):

Filtering statistics differ on the order of \( 10^{-6} \);
this is not pictured here.

Runge-Kutta step \( \Delta=10^{-2} \)

introduces small errors,
however, the residuals are unstructured in sign or magnitude.

Indicates random, unbiased numerical noise.

Shapiro–Wilk test has a p value of approximately 0.80.
Student’s t test with null hypothesis that the residuals have a mean of 0 has a p value of approximately 0.77.
We fail to reject the null hypothesis that the differences are distributed according to a mean zero Gaussian distribution.

Filter benchmark – Euler-Maruyama model-twin, fine discretization, Taylor truth-twin, fine discretization

Plot of of coarse Runge-Kutta filter simulation versus Taylor benchmark.

Generate ensemble with Euler-Maruyama:

step \( \Delta=10^{-3} \);
structure in the residuals;
artificial inflation in the forecast.

High diffusion:

RMSE and spread nearly identical to benchmark.

Impact of bias depends strongly on the model uncertainty.

Filter benchmark – Euler-Maruyama model-twin, coarse discretization, Taylor truth-twin, fine discretization

Plot of coarse Euler-Maruyama filter simulation versus Taylor benchmark.

Generate ensemble with Euler-Maruyama:

step \( \Delta=10^{-2} \);
structure in the residuals persists;
artificial inflation in the forecast.
Low diffusion: the error forces filter divergence.
High diffusion (\( s\geq 0.5 \)): difference with the benchmark is on the order of \( 10^{-2} \).

Filter benchmark – Runge-Kutta model-twin, coarse discretization, Taylor truth-twin, coarse discretization

Plot of of coarse Runge-Kutta filter simulation with coarse truth-twin versus Taylor benchmark.

Finally we test the benchmark configuration versus:

model-twin — Runge-Kutta step size \( \Delta=10^{-2} \); and
truth-twin — Taylor step size \( \Delta = 5 \times 10^{-3} \).

The expected discretization error is less than \( 10^{-3} \) over all diffusion regimes.
Inaccuracy in the truth-twin adds some bias in the residuals.
Student’s t test has p value on \( \mathcal{O}\left(10^{-4}\right) \), and we reject the hypothesis that the residuals are mean zero.
However, the practical difference in the statistics is negligible, with mean of the residuals at \( \mathcal{O}\left(10^{-4}\right) \).

Compared with:

model-twin — Euler-Maruyama; and
truth-twin — Taylor step size \( \Delta = 5 \times 10^{-3} \).
This relaxes the issues with Euler-Maruyama when the model-twin uses a step size of \( \Delta=10^{-3} \).
However, filter divergence still occurs when the model-twin uses a step size of \( \Delta=10^{-2} \).

Conclusions

We distinguish between strong and weak convergence, and its impact on the truth-twin and the model-twin respectively.
Euler-Maruyama:

introduces strong, systematic bias into twin experiments when the step size is greater than or equal to \( \Delta=10^{-2} \).
Effects depends strongly on model uncertainty in a twin experiment;
bias is sufficient to cause filter divergence in weak diffusion.

Runge-Kutta:

statistically robust solver.
\( \Delta \) on \( \mathcal{O}\left(10^{-3}\right) \), — virtually no difference with Taylor scheme.
Step \( \Delta=10^{-2} \) — introduces additional discretization error, but the error doesn’t strongly influence the RMSE or spread of the EnKF.

Conclusions – continued

We demonstrate in our work that an overall discretization error can be bounded by \( 10^{-3} \), using:

4-stage Runge-Kutta with step \( \Delta=10^{-2} \) — model-twin; and
order 2.0 Taylor with step \( \Delta=5\times 10^{-3} \) — truth-twin.

This forms a practical compromise, which our diagnostics demonstrate does not practically bias the outcomes of the filtering statistics.
We provide a computationally efficient framework for statistically robust twin experiments in the L96-s model.
We provide a novel derivation of ther order 2.0 Taylor scheme for the L96-s model.

This derivation has not previously appeared in the literature to the authors' knowledge, and we moreover provide benchmarks on the convergence of this and other schemes.

As an open question:

Can using weak integration schemes (which may not converge to any sample path) reduce the computational burden of ensemble-based forecasts in geophysical models?

If the goal of the forecast is to converge in distribution alone, this may be a viable alternative to traditional geophysical modeling paradigms.
Potential advantage:

Numerical precision can be targeted in the weak sense alone.